# NOTES

# Load package foreign 
# Load package MCMC (installed from .zip or make pkg-MCMCPack)

# Read data

data<-read.dta("jop_merge24.dta")
attach(data)
nvote<-(m+f)-v
x<-matrix(,nrow(data),2)
output<-matrix(,1,2)
XVAR<-cbind(urban, comp24, index, sumout, yss,constant)
XVAR<-XVAR*sqrt(n)

posterior<-MCMChierEXC(f,m,v, nvote, XVAR, burnin=2500, mcmc=7500, thin=3,
 delta0=2.0, delta1=2.0, nu1=1.0, nu0=1.0, verbose=TRUE)

for (i in 1:nrow(data)) {
x[i,1]<-mean(posterior[,i])
x[i,2]<-mean(posterior[,i+nrow(data)])
}

output<-matrix(,1,2)
output[1,1]<-weighted.mean(x[,1],data$n)
output[1,2]<-weighted.mean(x[,2],data$n)

 cat(" Pooled 1924 Data", "\n")
 cat(" MCMC iterations =", nrow(posterior), "\n")
# cat(" female turnout (weighted)=" , output[,1], "\n")
# cat(" male turnout (weighted)=", output[,2], "\n")
# cat("\n")

save.image("posterior_1924_jop.RData")

data$p0<-x[,1]
data$p1<-x[,2]
cat("RESULTS 1924", "\n" , "\n")
cat("Posterior median, 1924", "\n")
cat("Women: ", weighted.mean(data$p0,data$n), " Men: ", weighted.mean(data$p1,data$n))
cat("  Difference: ", (weighted.mean(data$p1,data$n)-weighted.mean(data$p0,data$n)), "\n")

cat("\n")

cat("Posterior medians, by state (1924)", "\n")
test <- subset(data, data$state == "1") 
cat("Massachusetts", "\n", "Women: ", weighted.mean(test$p0,test$n), " Men: ", weighted.mean(test$p1,test$n))
cat("  Difference: ", (weighted.mean(test$p1,test$n)-weighted.mean(test$p0,test$n)), "\n")
test <- subset(data, data$state == "3") 
cat("Connecticut ", "\n", "Women: ", weighted.mean(test$p0,test$n), " Men: ", weighted.mean(test$p1,test$n))
cat("  Difference: ", (weighted.mean(test$p1,test$n)-weighted.mean(test$p0,test$n)), "\n")
test <- subset(data, data$state == "13") 
cat("New York", "\n", "Women: ", weighted.mean(test$p0,test$n), " Men: ", weighted.mean(test$p1,test$n))
cat("  Difference: ", (weighted.mean(test$p1,test$n)-weighted.mean(test$p0,test$n)), "\n")
test <- subset(data, data$state == "21") 
cat("Illinois ", "\n", "Women: ", weighted.mean(test$p0,test$n), " Men: ", weighted.mean(test$p1,test$n))
cat("  Difference: ", (weighted.mean(test$p1,test$n)-weighted.mean(test$p0,test$n)), "\n")
test <- subset(data, data$state == "23") 
cat("Michigan", "\n", "Women: ", weighted.mean(test$p0,test$n), " Men: ", weighted.mean(test$p1,test$n))
cat("  Difference: ", (weighted.mean(test$p1,test$n)-weighted.mean(test$p0,test$n)), "\n","\n")

cat("Posterior medians, by urban (1924)", "\n")
test <- subset(data, data$u == "0") 
cat("Rural", "\n", "Women: ", weighted.mean(test$p0,test$n), " Men: ", weighted.mean(test$p1,test$n))
cat("  Difference: ", (weighted.mean(test$p1,test$n)-weighted.mean(test$p0,test$n)), "\n")
test <- subset(data, data$u == "1") 
cat("Urban", "\n", "Women: ", weighted.mean(test$p0,test$n), " Men: ", weighted.mean(test$p1,test$n))
cat("  Difference: ", (weighted.mean(test$p1,test$n)-weighted.mean(test$p0,test$n)), "\n", "\n")

cat("Regression coefficients", "\n")

b10<-median(posterior[,1+(nrow(data)*2)])
b20<-median(posterior[,2+(nrow(data)*2)])
b30<-median(posterior[,3+(nrow(data)*2)])
b40<-median(posterior[,4+(nrow(data)*2)])
b50<-median(posterior[,5+(nrow(data)*2)])
b00<-median(posterior[,6+(nrow(data)*2)])

b11<-median(posterior[,7+(nrow(data)*2)])
b21<-median(posterior[,8+(nrow(data)*2)])
b31<-median(posterior[,9+(nrow(data)*2)])
b41<-median(posterior[,10+(nrow(data)*2)])
b51<-median(posterior[,11+(nrow(data)*2)])
b01<-median(posterior[,12+(nrow(data)*2)])

# COMPUTE REGION OF HIGHEST POSTERIOR DENSITY
hpd1<-matrix(nrow=6,ncol=2)
for (j in (nrow(data)*2+1):(nrow(data)*2+6)){
z<-matrix(nrow=500, ncol=2)
for (i in 1:500) {
z[i,1]<-i
z[i,2]<-quantile(posterior[,j],probs=c(0.95+((0.05*(i/500)))))-quantile(posterior[,j], probs=c(0.0+((0.05*(i/500)))))
}
# Identity quantile pair with smallest segment
a<-z[,1][z[,2]==min(z[,2])]
# compute quantiles for that setting of i
hpd1[j-(nrow(data)*2),]<-quantile(posterior[,j],probs=c(0.00+((0.05*a/500)),0.95+(0.05*a/500)))
}

# COMPUTE REGION OF HIGHEST POSTERIOR DENSITY
hpd2<-matrix(nrow=6,ncol=2)
for (j in (nrow(data)*2+7):(nrow(data)*2+12)){
z<-matrix(nrow=500, ncol=2)
for (i in 1:500) {
z[i,1]<-i
z[i,2]<-quantile(posterior[,j],probs=c(0.95+((0.05*(i/500)))))-quantile(posterior[,j], probs=c(0.0+((0.05*(i/500)))))
}
# Identity quantile pair with smallest segment
a<-z[,1][z[,2]==min(z[,2])]
# compute quantiles for that setting of i
hpd2[j-(nrow(data)*2)-6,]<-quantile(posterior[,j],probs=c(0.00+((0.05*a/500)),0.95+(0.05*a/500)))
}

cat("FEMALE, 1924", "\n")
cat("Urban               ", b10, "  ", hpd1[1,],"\n")
cat("Comp 24             ", b20, "  ", hpd1[2,],"\n")
cat("Electoral Laws      ", b30, "  ", hpd1[3,],"\n")
cat("Outside activity    ", b40, "  ", hpd1[4,],"\n")
cat("Years Since Suffrage", b50, "  ", hpd1[5,],"\n")
cat("Constant,           ", b00, "  ", hpd1[6,],"\n")

cat("MALE, 1924", "\n")
cat("Urban               ", b11, "  ", hpd2[1,],"\n")
cat("Comp 24            ", b21, "  ", hpd2[2,],"\n")
cat("Electoral Laws      ", b31, "  ", hpd2[3,],"\n")
cat("Outside activity    ", b41, "  ", hpd2[4,],"\n")
cat("Years Since Suffrage", b51, "  ", hpd2[5,],"\n")
cat("Constant            ", b01, "  ", hpd2[6,],"\n", "\n")

cat("Effects of independent variables", "\n")
attach(data)
pu0<-(1/(1+exp(-(b00+b10*min(urban)+b20*mean(comp)+b30*mean(index)+b40*mean(sumout)+b50*mean(yss)))))
pu1<-(1/(1+exp(-(b00+b10*max(urban)+b20*mean(comp)+b30*mean(index)+b40*mean(sumout)+b50*mean(yss)))))
pc0<-(1/(1+exp(-(b00+b10*mean(urban)+b20*min(comp)+b30*mean(index)+b40*mean(sumout)+b50*mean(yss)))))
pc1<-(1/(1+exp(-(b00+b10*mean(urban)+b20*max(comp)+b30*mean(index)+b40*mean(sumout)+b50*mean(yss)))))
pe0<-(1/(1+exp(-(b00+b10*mean(urban)+b20*mean(comp)+b30*min(index)+b40*mean(sumout)+b50*mean(yss)))))
pe1<-(1/(1+exp(-(b00+b10*mean(urban)+b20*mean(comp)+b30*max(index)+b40*mean(sumout)+b50*mean(yss)))))
po0<-(1/(1+exp(-(b00+b10*mean(urban)+b20*mean(comp)+b30*mean(index)+b40*min(sumout)+b50*mean(yss)))))
po1<-(1/(1+exp(-(b00+b10*mean(urban)+b20*mean(comp)+b30*mean(index)+b40*max(sumout)+b50*mean(yss)))))
py0<-(1/(1+exp(-(b00+b10*mean(urban)+b20*mean(comp)+b30*mean(index)+b40*mean(sumout)+b50*min(yss)))))
py1<-(1/(1+exp(-(b00+b10*mean(urban)+b20*mean(comp)+b30*mean(index)+b40*mean(sumout)+b50*max(yss)))))
cat("FEMALE, 1924", "\n")
cat("Urban               ", pu1-pu0, "\n")
cat("Comp 24             ", pc1-pc0, "\n")
cat("Electoral Laws      ", pe1-pe0, "\n")
cat("Outside activity    ", po1-po0, "\n")
cat("Years Since Suffrage", py1-py0, "\n","\n")



pu0<-(1/(1+exp(-(b01+b11*min(urban)+b21*mean(comp)+b31*mean(index)+b41*mean(sumout)+b51*mean(yss)))))
pu1<-(1/(1+exp(-(b01+b11*max(urban)+b21*mean(comp)+b31*mean(index)+b41*mean(sumout)+b51*mean(yss)))))
pc0<-(1/(1+exp(-(b01+b11*mean(urban)+b21*min(comp)+b31*mean(index)+b41*mean(sumout)+b51*mean(yss)))))
pc1<-(1/(1+exp(-(b01+b11*mean(urban)+b21*max(comp)+b31*mean(index)+b41*mean(sumout)+b51*mean(yss)))))
pe0<-(1/(1+exp(-(b01+b11*mean(urban)+b21*mean(comp)+b31*min(index)+b41*mean(sumout)+b51*mean(yss)))))
pe1<-(1/(1+exp(-(b01+b11*mean(urban)+b21*mean(comp)+b31*max(index)+b41*mean(sumout)+b51*mean(yss)))))
po0<-(1/(1+exp(-(b01+b11*mean(urban)+b21*mean(comp)+b31*mean(index)+b41*min(sumout)+b51*mean(yss)))))
po1<-(1/(1+exp(-(b01+b11*mean(urban)+b21*mean(comp)+b31*mean(index)+b41*max(sumout)+b51*mean(yss)))))
py0<-(1/(1+exp(-(b01+b11*mean(urban)+b21*mean(comp)+b31*mean(index)+b41*mean(sumout)+b51*min(yss)))))
py1<-(1/(1+exp(-(b01+b11*mean(urban)+b21*mean(comp)+b31*mean(index)+b41*mean(sumout)+b51*max(yss)))))
cat("MALE, 1924", "\n")
cat("Urban               ", pu1-pu0, "\n")
cat("Comp 24            ", pc1-pc0, "\n")
cat("Electoral Laws      ", pe1-pe0, "\n")
cat("Outside activity    ", po1-po0, "\n")
cat("Years Since Suffrage", py1-py0, "\n","\n")


















detach(data)

